library(topGO)
library(stats)
library(knitr)


setwd('/Users/mmayers/projects/n15_mice/data/')

data.dir <- '../data'

gene2GOID <- readMappings('clusterID2GO.map')
unenr_pvals <- read.csv(file.path(data.dir, 'unenriched_pvals.csv'))
unenr_pvals$log2FoldChange <- log2(unenr_pvals$ratio)
unenr_pvals$padj <- p.adjust(unenr_pvals$p_value, method="BH")
GO2geneMF <- annFUN.gene2GO('MF', gene2GO=gene2GOID)
GO2geneBP <- annFUN.gene2GO('BP', gene2GO=gene2GOID)

do some stuff with fold changes

unenr_pvals$adjFC <- unenr_pvals$log2FoldChange
n <- sum(unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2)
unenr_pvals[unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2, ]$adjFC <- rep(0, n)
#set.seed(2222)
#unenr_pvals[unenr_pvals$padj > .05 & abs(unenr_pvals$log2FoldChange) > 2, ]$adjFC <- runif(n, -1, 1)
hist(unenr_pvals$log2FoldChange, 15)

hist(unenr_pvals$adjFC, 15)

Metadata on Cluster Info - Use this for filtering of certain types of proteins

sampleTable <- read.csv(file.path(data.dir, 'filt_metadata.csv'), row.names=1)
sampleTable$technical <- as.logical(sampleTable$technical)
locusTable <- read.csv(file.path(data.dir, 'loci_annot.csv'))

Look at Bacterial Proteins only

Mouse Human subset seems to be a bust, can’t figure out good way to account for tiny background)

tmp <- locusTable$mouse_human == 'False' & locusTable$lca != '35823' # Filtering out anthrospira differences 35823
tmp[is.na(tmp)] <- TRUE
bac.prots <- locusTable[tmp, ]$X
bac.pvals <- unenr_pvals[unenr_pvals$X %in% bac.prots, ]
geneList.bac <- bac.pvals$padj
geneList.bac <- bac.pvals$adjFC
names(geneList.bac) <- bac.pvals$X

topDiffGenes <- function(score) {return(score >= 2)} # scoring function, anything with > 2 logfc is sig anyway

TCell vs RAG - Unenriched

geneList.bac <- -1*geneList.bac

MF first

ontology <- 'MF'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0004619 phosphoglycerate mutase activity 42 0 4.77 1.4e-26
GO:0000287 magnesium ion binding 46 0 5.23 8.6e-09
GO:0005506 iron ion binding 80 0 9.09 7.2e-07
GO:0016491 oxidoreductase activity 291 15 33.07 1.3e-06
GO:0008184 glycogen phosphorylase activity 14 0 1.59 1.6e-06
GO:0004634 phosphopyruvate hydratase activity 34 0 3.86 5.0e-06
GO:0046933 proton-transporting ATP synthase activit… 14 1 1.59 4.1e-05
GO:0016639 oxidoreductase activity, acting on the C… 29 0 3.30 4.3e-05
GO:0004618 phosphoglycerate kinase activity 27 3 3.07 6.3e-05
GO:0003677 DNA binding 61 0 6.93 6.8e-05
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}
## [1] "phosphoglycerate mutase activity"

## [1] "magnesium ion binding"

## [1] "iron ion binding"

## [1] "oxidoreductase activity"

## [1] "glycogen phosphorylase activity"

## [1] "phosphopyruvate hydratase activity"

## [1] "proton-transporting ATP synthase activit..."

## [1] "oxidoreductase activity, acting on the C..."

## [1] "phosphoglycerate kinase activity"

## [1] "DNA binding"

#showGroupDensity(GOdata.bac, "GO:0008184", ranks=TRUE) #Glycogen phosphorlase
#showGroupDensity(GOdata.bac, "GO:0004618", ranks=TRUE) #phosphoglycerate kinase
ontology <- 'BP'

BP - Tcell vs RAG: Unenriched

GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0006096 glycolytic process 148 25 15.31 2.3e-18
GO:0005975 carbohydrate metabolic process 320 48 33.11 4.0e-14
GO:0008152 metabolic process 1017 92 105.24 2.1e-06
GO:0055114 oxidation-reduction process 312 28 32.29 1.8e-05
GO:0015986 ATP synthesis coupled proton transport 15 1 1.55 0.00028
GO:0006520 cellular amino acid metabolic process 78 0 8.07 0.00031
GO:0006108 malate metabolic process 6 0 0.62 0.00037
GO:0008654 phospholipid biosynthetic process 13 0 1.35 0.00072
GO:0006021 inositol biosynthetic process 13 0 1.35 0.00072
GO:0006351 transcription, DNA-templated 63 0 6.52 0.00077
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}
## [1] "glycolytic process"

## [1] "carbohydrate metabolic process"

## [1] "metabolic process"

## [1] "oxidation-reduction process"

## [1] "ATP synthesis coupled proton transport"

## [1] "cellular amino acid metabolic process"

## [1] "malate metabolic process"

## [1] "phospholipid biosynthetic process"

## [1] "inositol biosynthetic process"

## [1] "transcription, DNA-templated"

#showGroupDensity(GOdata.bac, "GO:0006096", ranks=TRUE) # Glycolitc Process
#showGroupDensity(GOdata.bac, "GO:0005975", ranks=TRUE) # Carbohydydrate metab
# Mostly Carb metabolism?  lipid biosynthesis?Nah, those look crappy

RAG vs Tcell - Unenriched

geneList.bac <- -1*geneList.bac

MF

ontology <- 'MF'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0005524 ATP binding 233 2 4.38 1.7e-24
GO:0003872 6-phosphofructokinase activity 25 0 0.47 1.8e-14
GO:0004478 methionine adenosyltransferase activity 29 0 0.55 8.1e-13
GO:0005198 structural molecule activity 151 0 2.84 5.5e-12
GO:0050242 pyruvate, phosphate dikinase activity 74 0 1.39 6.9e-12
GO:0016887 ATPase activity 34 0 0.64 2.5e-11
GO:0008878 glucose-1-phosphate adenylyltransferase … 15 0 0.28 7.0e-11
GO:0016301 kinase activity 128 1 2.41 5.3e-09
GO:0003924 GTPase activity 77 1 1.45 4.1e-08
GO:0003824 catalytic activity 1023 22 19.24 1.7e-06
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}
## [1] "ATP binding"

## [1] "6-phosphofructokinase activity"

## [1] "methionine adenosyltransferase activity"

## [1] "structural molecule activity"

## [1] "pyruvate, phosphate dikinase activity"

## [1] "ATPase activity"

## [1] "glucose-1-phosphate adenylyltransferase ..."

## [1] "kinase activity"

## [1] "GTPase activity"

## [1] "catalytic activity"

BP - RAG vs Tcell: Unenriched

ontology <- 'BP'
GOdata.bac <- new('topGOdata', ontology = ontology, allGenes = geneList.bac, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.ks <- runTest(GOdata.bac, statistic = 'ks')
a <- GenTable(GOdata.bac, result.bac.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0006556 S-adenosylmethionine biosynthetic proces… 29 0 0.48 6.8e-14
GO:0071973 bacterial-type flagellum-dependent cell … 95 0 1.57 3.8e-13
GO:0006002 fructose 6-phosphate metabolic process 22 0 0.36 7.4e-13
GO:0006090 pyruvate metabolic process 222 0 3.67 1.1e-12
GO:0016310 phosphorylation 229 1 3.78 1.2e-12
GO:0005978 glycogen biosynthetic process 15 0 0.25 4.1e-11
GO:0009401 phosphoenolpyruvate-dependent sugar phos… 13 0 0.21 4.6e-07
GO:0044267 cellular protein metabolic process 130 2 2.15 1.1e-06
GO:0006810 transport 61 0 1.01 4.0e-06
GO:0055114 oxidation-reduction process 312 4 5.15 4.6e-05
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac, term, ranks=TRUE))
}
## [1] "S-adenosylmethionine biosynthetic proces..."

## [1] "bacterial-type flagellum-dependent cell ..."

## [1] "fructose 6-phosphate metabolic process"

## [1] "pyruvate metabolic process"

## [1] "phosphorylation"

## [1] "glycogen biosynthetic process"

## [1] "phosphoenolpyruvate-dependent sugar phos..."

## [1] "cellular protein metabolic process"

## [1] "transport"

## [1] "oxidation-reduction process"

BioGlyCMK Rag vs Tcell

enr_pvals <- read.csv(file.path(data.dir, 'enriched_pvals.csv'))
enr_pvals$log2FoldChange <- log2(enr_pvals$ratio)
enr_pvals$padj <- p.adjust(enr_pvals$p_value, method="BH")

Fix fold change of insignifcant proteins

enr_pvals$adjFC <- enr_pvals$log2FoldChange
n <- sum(enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2)
#set.seed(1234)
#enr_pvals[enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2, ]$adjFC <- runif(n, -1, 1)
enr_pvals[enr_pvals$padj > .05 & abs(enr_pvals$log2FoldChange) > 2, ]$adjFC <- rep(0, n)
hist(enr_pvals$log2FoldChange, 15)

hist(enr_pvals$adjFC, 15)

Bacterial proteins only

bac.enr.pvals <- enr_pvals[enr_pvals$X %in% bac.prots, ]
#geneList.bac <- bac.pvals$padj
geneList.bac.enr <- bac.enr.pvals$adjFC
names(geneList.bac.enr) <- bac.enr.pvals$X

Tcell vs RAG - Enriched

geneList.bac.enr <- -1*geneList.bac.enr # KS goes from - to +
                                        # do this to put upreg in Tcell first

MF

ontology <- 'MF'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0016639 oxidoreductase activity, acting on the C… 70 0 0.45 2.5e-07
GO:0016491 oxidoreductase activity 367 2 2.36 6.1e-07
GO:0016829 lyase activity 160 1 1.03 1.2e-06
GO:0051287 NAD binding 22 0 0.14 5.1e-06
GO:0016462 pyrophosphatase activity 50 0 0.32 6.6e-06
GO:0050661 NADP binding 20 0 0.13 9.4e-06
GO:0005515 protein binding 16 0 0.10 0.00020
GO:0005215 transporter activity 25 0 0.16 0.00025
GO:0004197 cysteine-type endopeptidase activity 7 0 0.04 0.00135
GO:0005506 iron ion binding 52 0 0.33 0.00155
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}
## [1] "oxidoreductase activity, acting on the C..."

## [1] "oxidoreductase activity"

## [1] "lyase activity"

## [1] "NAD binding"

## [1] "pyrophosphatase activity"

## [1] "NADP binding"

## [1] "protein binding"

## [1] "transporter activity"

## [1] "cysteine-type endopeptidase activity"

## [1] "iron ion binding"

BP Tcell vs RAG: BioGlyCMK

ontology <- 'BP'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0006520 cellular amino acid metabolic process 140 0 0.65 7.7e-13
GO:0055114 oxidation-reduction process 369 2 1.72 7.6e-10
GO:0006508 proteolysis 19 0 0.09 7.4e-07
GO:0006814 sodium ion transport 18 0 0.08 2.0e-06
GO:0006006 glucose metabolic process 53 0 0.25 6.4e-06
GO:0044267 cellular protein metabolic process 35 0 0.16 2.5e-05
GO:0042026 protein refolding 19 0 0.09 0.00067
GO:0005980 glycogen catabolic process 4 0 0.02 0.00150
GO:0006555 methionine metabolic process 4 0 0.02 0.00462
GO:0006177 GMP biosynthetic process 3 0 0.01 0.00569
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}
## [1] "cellular amino acid metabolic process"

## [1] "oxidation-reduction process"

## [1] "proteolysis"

## [1] "sodium ion transport"

## [1] "glucose metabolic process"

## [1] "cellular protein metabolic process"

## [1] "protein refolding"

## [1] "glycogen catabolic process"

## [1] "methionine metabolic process"

## [1] "GMP biosynthetic process"

#showGroupDensity(GOdata.bac.enr, "GO:0006508", ranks = TRUE) #proteolysis

RAG vs TCell - BioGlyCMK

ontology <- 'MF'

MF

geneList.bac.enr <- -1*geneList.bac.enr 

GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0050242 pyruvate, phosphate dikinase activity 55 0 0.65 1.3e-25
GO:0016301 kinase activity 61 0 0.72 3.2e-19
GO:0016903 oxidoreductase activity, acting on the a… 141 0 1.66 9.1e-13
GO:0051536 iron-sulfur cluster binding 45 0 0.53 2.7e-11
GO:0016887 ATPase activity 33 0 0.39 9.0e-11
GO:0005524 ATP binding 149 0 1.75 4.7e-10
GO:0004634 phosphopyruvate hydratase activity 44 0 0.52 2.2e-08
GO:0016820 hydrolase activity, acting on acid anhyd… 16 0 0.19 6.8e-07
GO:0000287 magnesium ion binding 57 1 0.67 1.0e-05
GO:0016747 transferase activity, transferring acyl … 26 0 0.31 1.6e-05
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}
## [1] "pyruvate, phosphate dikinase activity"

## [1] "kinase activity"

## [1] "oxidoreductase activity, acting on the a..."

## [1] "iron-sulfur cluster binding"

## [1] "ATPase activity"

## [1] "ATP binding"

## [1] "phosphopyruvate hydratase activity"

## [1] "hydrolase activity, acting on acid anhyd..."

## [1] "magnesium ion binding"

## [1] "transferase activity, transferring acyl ..."

BP - RAG vs Tcell: BioGlyCMK

ontology <- 'BP'
GOdata.bac.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.bac.enr, geneSel = topDiffGenes, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.bac.enr.ks <- runTest(GOdata.bac.enr, statistic = 'ks')
a <- GenTable(GOdata.bac.enr, result.bac.enr.ks, topNodes = 10)
kable(a)
GO.ID Term Annotated Significant Expected result1
GO:0006090 pyruvate metabolic process 132 4 1.54 1.2e-25
GO:0016310 phosphorylation 133 4 1.55 4.4e-25
GO:0022900 electron transport chain 16 0 0.19 2.9e-05
GO:0005978 glycogen biosynthetic process 4 0 0.05 0.00044
GO:0006096 glycolytic process 77 4 0.90 0.00046
GO:0016226 iron-sulfur cluster assembly 12 0 0.14 0.00054
GO:0006730 one-carbon metabolic process 7 0 0.08 0.00286
GO:0006810 transport 51 0 0.60 0.00352
GO:0005975 carbohydrate metabolic process 224 4 2.61 0.00414
GO:0006091 generation of precursor metabolites and … 116 4 1.35 0.00616
for (term in a$GO.ID){
  print(a[a$GO.ID == term, ]$Term)
  print(showGroupDensity(GOdata.bac.enr, term, ranks=TRUE))
}
## [1] "pyruvate metabolic process"

## [1] "phosphorylation"

## [1] "electron transport chain"

## [1] "glycogen biosynthetic process"

## [1] "glycolytic process"

## [1] "iron-sulfur cluster assembly"

## [1] "one-carbon metabolic process"

## [1] "transport"

## [1] "carbohydrate metabolic process"

## [1] "generation of precursor metabolites and ..."

#showGroupDensity(GOdata.bac.enr, "GO:0006091", ranks=TRUE) #generation of precursor metabolites
#showGroupDensity(GOdata.bac.enr, "GO:0006090", ranks=TRUE) # pyruvate metabolic process

Enriched vs Unenriched

groups <- read.csv(file.path(data.dir, 'groups.csv'))
groups.rt <- groups[groups$RT_Enriched == 'True' | groups$RT_Unenriched == 'True', ]
geneList.rt.enr <- groups.rt$RT_Enriched == 'True' & groups.rt$RT_Unenriched == 'False'
geneList.rt.enr <- as.numeric(geneList.rt.enr)
names(geneList.rt.enr) <- groups.rt$X

sigFunc <- function(score) {return(score > 0.5)}

MF - Enr vs Unenr: Tcell

ontology <- 'MF'

GOdata.rt.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.enr.fisher <- runTest(GOdata.rt.enr, statistic = 'fisher')
kable(GenTable(GOdata.rt.enr, result.rt.enr.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0016620 oxidoreductase activity, acting on the a… 146 76 30.70 1.2e-14
GO:0004816 asparagine-tRNA ligase activity 18 17 3.78 3.8e-11
GO:0003824 catalytic activity 2849 644 599.07 5.0e-08
GO:0008234 cysteine-type peptidase activity 23 19 4.84 3.4e-06
GO:0046872 metal ion binding 472 83 99.25 1.3e-05
GO:0008774 acetaldehyde dehydrogenase (acetylating)… 17 12 3.57 1.5e-05
GO:0004197 cysteine-type endopeptidase activity 15 11 3.15 2.0e-05
GO:0004022 alcohol dehydrogenase (NAD) activity 18 12 3.78 3.6e-05
GO:0004375 glycine dehydrogenase (decarboxylating) … 6 6 1.26 8.5e-05
GO:0004825 methionine-tRNA ligase activity 12 9 2.52 9.1e-05

BP - Enr vs Unenr: Tcell

ontology <- 'BP'

GOdata.rt.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.enr.fisher <- runTest(GOdata.rt.enr, statistic = 'fisher')
kable(GenTable(GOdata.rt.enr, result.rt.enr.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0006421 asparaginyl-tRNA aminoacylation 18 17 3.71 2.7e-11
GO:0008152 metabolic process 2834 603 583.98 1.0e-10
GO:0015976 carbon utilization 17 12 3.50 1.2e-05
GO:0006508 proteolysis 89 36 18.34 3.9e-05
GO:0006431 methionyl-tRNA aminoacylation 12 9 2.47 7.7e-05
GO:0006066 alcohol metabolic process 42 13 8.65 0.00014
GO:0009082 branched-chain amino acid biosynthetic p… 24 13 4.95 0.00028
GO:0009081 branched-chain amino acid metabolic proc… 29 18 5.98 0.00034
GO:0006099 tricarboxylic acid cycle 5 5 1.03 0.00037
GO:0006546 glycine catabolic process 8 7 1.65 0.00043

RAG Enr vs Unenr

groups.rag <- groups[groups$RAG_Enriched == 'True' | groups$RAG_Unenriched == 'True', ]
geneList.rag.enr <- groups.rag$RAG_Enriched == 'True' & groups.rag$RAG_Unenriched == 'False'
geneList.rag.enr <- as.numeric(geneList.rag.enr)
names(geneList.rag.enr) <- groups.rag$X

MF - Enr vs Unenr: RAG

ontology <- 'MF'

GOdata.rag.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.enr.fisher <- runTest(GOdata.rag.enr, statistic = 'fisher')
kable(GenTable(GOdata.rag.enr, result.rag.enr.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0016620 oxidoreductase activity, acting on the a… 142 77 29.01 4.7e-17
GO:0004816 asparagine-tRNA ligase activity 18 17 3.68 2.3e-11
GO:0003824 catalytic activity 2697 595 550.94 1.2e-08
GO:0008774 acetaldehyde dehydrogenase (acetylating)… 17 12 3.47 1.1e-05
GO:0008234 cysteine-type peptidase activity 22 18 4.49 1.3e-05
GO:0004197 cysteine-type endopeptidase activity 15 11 3.06 1.5e-05
GO:0004022 alcohol dehydrogenase (NAD) activity 18 12 3.68 2.6e-05
GO:0004375 glycine dehydrogenase (decarboxylating) … 6 6 1.23 7.1e-05
GO:0004825 methionine-tRNA ligase activity 12 9 2.45 7.2e-05
GO:0004019 adenylosuccinate synthase activity 8 7 1.63 9.5e-05

BP - Enr vs Unenr: RAG

ontology <- 'BP'

GOdata.rag.enr <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.enr, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.enr.fisher <- runTest(GOdata.rag.enr, statistic = 'fisher')
kable(GenTable(GOdata.rag.enr, result.rag.enr.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0006421 asparaginyl-tRNA aminoacylation 18 17 3.55 1.2e-11
GO:0008152 metabolic process 2674 547 527.01 3.8e-10
GO:0015976 carbon utilization 17 12 3.35 7.3e-06
GO:0006508 proteolysis 84 34 16.56 2.7e-05
GO:0006431 methionyl-tRNA aminoacylation 12 9 2.37 5.3e-05
GO:0055114 oxidation-reduction process 867 207 170.88 6.4e-05
GO:0006066 alcohol metabolic process 42 13 8.28 8.7e-05
GO:0009082 branched-chain amino acid biosynthetic p… 24 13 4.73 0.00017
GO:0009081 branched-chain amino acid metabolic proc… 28 17 5.52 0.00141
GO:0006546 glycine catabolic process 7 6 1.38 0.00146

Unenriched vs Enriched

geneList.rt.un <- groups.rt$RT_Enriched == 'False' & groups.rt$RT_Unenriched == 'True'
geneList.rt.un <- as.numeric(geneList.rt.un)
names(geneList.rt.un) <- groups.rt$X

MF - Unenr vs Enr: Tcell

ontology <- 'MF'

GOdata.rt.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.un.fisher <- runTest(GOdata.rt.un, statistic = 'fisher')
kable(GenTable(GOdata.rt.un, result.rt.un.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0003735 structural constituent of ribosome 153 132 47.98 < 1e-30
GO:0005198 structural molecule activity 298 246 93.44 < 1e-30
GO:0003677 DNA binding 137 111 42.96 < 1e-30
GO:0003723 RNA binding 216 145 67.73 < 1e-30
GO:0032549 ribonucleoside binding 861 279 269.98 < 1e-30
GO:0003899 DNA-directed RNA polymerase activity 137 107 42.96 1.7e-30
GO:0004619 phosphoglycerate mutase activity 44 40 13.80 1.4e-16
GO:0046933 proton-transporting ATP synthase activit… 32 30 10.03 1.4e-13
GO:0016740 transferase activity 869 284 272.49 0.00019
GO:0003924 GTPase activity 119 55 37.31 0.00039

BP - Unenr vs Enr: Tcell

ontology <- 'BP'

GOdata.rt.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rt.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rt.un.fisher <- runTest(GOdata.rt.un, statistic = 'fisher')
kable(GenTable(GOdata.rt.un, result.rt.un.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0071973 bacterial-type flagellum-dependent cell … 112 112 36.01 < 1e-30
GO:0006412 translation 261 173 83.92 < 1e-30
GO:0006351 transcription, DNA-templated 152 116 48.87 9.1e-30
GO:0015986 ATP synthesis coupled proton transport 34 30 10.93 1.3e-11
GO:0000162 tryptophan biosynthetic process 5 5 1.61 0.0034
GO:0006438 valyl-tRNA aminoacylation 5 5 1.61 0.0034
GO:0006556 S-adenosylmethionine biosynthetic proces… 84 39 27.01 0.0039
GO:0015991 ATP hydrolysis coupled proton transport 37 20 11.90 0.0045
GO:0006801 superoxide metabolic process 7 6 2.25 0.0056
GO:0006563 L-serine metabolic process 15 11 4.82 0.0066

RAG Unenriched vs Enriched

geneList.rag.un <- groups.rag$RAG_Enriched == 'False' & groups.rag$RAG_Unenriched == 'True'
geneList.rag.un <- as.numeric(geneList.rag.un)
names(geneList.rag.un) <- groups.rag$X

MF - Unenr vs Enr: RAG

ontology <- 'MF'

GOdata.rag.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.un.fisher <- runTest(GOdata.rag.un, statistic = 'fisher')
kable(GenTable(GOdata.rag.un, result.rag.un.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0005198 structural molecule activity 271 244 92.18 < 1e-30
GO:0003735 structural constituent of ribosome 139 124 47.28 < 1e-30
GO:0003677 DNA binding 136 111 46.26 < 1e-30
GO:0032549 ribonucleoside binding 835 285 284.04 2.2e-30
GO:0003723 RNA binding 202 146 68.71 8.6e-29
GO:0003899 DNA-directed RNA polymerase activity 136 107 46.26 1.3e-27
GO:0004619 phosphoglycerate mutase activity 43 42 14.63 3.7e-19
GO:0046933 proton-transporting ATP synthase activit… 33 31 11.23 5.6e-13
GO:0016661 oxidoreductase activity, acting on other… 23 22 7.82 6.8e-10
GO:0003924 GTPase activity 113 62 38.44 3.0e-06

BP - Unenr vs Enr: RAG

ontology <- 'BP'

GOdata.rag.un <- new('topGOdata', ontology = ontology, allGenes = geneList.rag.un, geneSel = sigFunc, annot = annFUN.gene2GO, gene2GO = gene2GOID)
result.rag.un.fisher <- runTest(GOdata.rag.un, statistic = 'fisher')
kable(GenTable(GOdata.rag.un, result.rag.un.fisher, topNodes = 10))
GO.ID Term Annotated Significant Expected result1
GO:0071973 bacterial-type flagellum-dependent cell … 112 112 38.53 < 1e-30
GO:0006412 translation 237 167 81.54 < 1e-30
GO:0006351 transcription, DNA-templated 149 114 51.26 2.2e-27
GO:0015986 ATP synthesis coupled proton transport 35 31 12.04 3.4e-11
GO:0015991 ATP hydrolysis coupled proton transport 32 21 11.01 0.00029
GO:0006414 translational elongation 56 31 19.27 0.00094
GO:0000162 tryptophan biosynthetic process 5 5 1.72 0.00479
GO:0006801 superoxide metabolic process 7 6 2.41 0.00812
GO:0006563 L-serine metabolic process 15 11 5.16 0.01086
GO:0006556 S-adenosylmethionine biosynthetic proces… 84 39 28.90 0.01385